home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / DE118I.ZIP / PRECESS.C < prev    next >
C/C++ Source or Header  |  1993-01-30  |  13KB  |  473 lines
  1. /* Precession of the equinox and ecliptic
  2.  * from epoch Julian date J to or from J2000.0
  3.  *
  4.  * Program by Steve Moshier.
  5.  *
  6.  *
  7.  * The IAU formula is not used in this version.
  8.  * IAU Coefficients are from:
  9.  * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
  10.  * "Expressions for the Precession Quantities Based upon the IAU
  11.  * (1976) System of Astronomical Constants,"  Astronomy and
  12.  * Astrophysics 58, 1-16 (1977).
  13.  *
  14.  * Newer formulas that cover a much longer time span are from:
  15.  * J. Laskar, "Secular terms of classical planetary theories
  16.  * using the results of general theory," Astronomy and Astrophysics
  17.  * 157, 59070 (1986).
  18.  *
  19.  * See also:
  20.  * P. Bretagnon and G. Francou, "Planetary theories in rectangular
  21.  * and spherical variables. VSOP87 solutions," Astronomy and
  22.  * Astrophysics 202, 309-315 (1988).
  23.  *
  24.  * Laskar's expansions are said by Bretagnon and Francou
  25.  * to have "a precision of about 1" over 10000 years before
  26.  * and after J2000.0 in so far as the precession constants p^0_A
  27.  * and epsilon^0_A are perfectly known."
  28.  *
  29.  * Bretagnon and Francou's expansions for the node and inclination
  30.  * of the ecliptic were derived from Laskar's data but were truncated
  31.  * after the term in T**6. I have recomputed these expansions from
  32.  * Laskar's data, retaining powers up to T**10 in the result.
  33.  *
  34.  * The following table indicates the differences between the result
  35.  * of the IAU formula and Laskar's formula using four different test
  36.  * vectors, checking at J2000 plus and minus the indicated number
  37.  * of years.
  38.  *
  39.  *   Years       Arc
  40.  * from J2000  Seconds
  41.  * ----------  -------
  42.  *        0      0
  43.  *      100    .006    
  44.  *      200     .006
  45.  *      500     .015
  46.  *     1000     .28
  47.  *     2000    6.4
  48.  *     3000   38.
  49.  *    10000 9400.
  50.  */
  51.  
  52. #ifndef NOINCS
  53. #include "mconf.h"
  54. #include "prec.h"
  55. #endif
  56. #if LDOUBLE
  57. #if UNK
  58. /* Precession coefficients taken from Laskar's paper: */
  59. static DOUBLE pAcof[10] = {
  60. -8.66000000000000000000E-10L,
  61. -4.75900000000000000000E-8L,
  62.  2.42400000000000000000E-7L,
  63.  1.30950000000000000000E-5L,
  64.  1.74510000000000000000E-4L,
  65. -1.80550000000000000000E-3L,
  66. -2.35316000000000000000E-1L,
  67.  7.73200000000000000000E-2L,
  68.  1.11197100000000000000E2L,
  69.  5.02909660000000000000E4L
  70. };
  71. /* Node and inclination of the earth's orbit computed from
  72.  * Laskar's data as done in Bretagnon and Francou's paper:
  73.  */
  74. static DOUBLE nodecof[11] = {
  75.  6.64020000000000000000E-16L,
  76. -2.69151000000000000000E-15L,
  77. -1.54702100000000000000E-12L,
  78.  7.52131300000000000000E-12L,
  79.  6.31901310000000000000E-10L,
  80. -3.48388152000000000000E-9L,
  81. -1.81306589600000000000E-7L,
  82.  2.75036225000000000000E-8L,
  83.  7.43945314260000000000E-5L,
  84. -4.20786043170000000000E-2L,
  85.  3.05211265497500000000E0L
  86. };
  87. static DOUBLE inclcof[11] = {
  88.  1.21470000000000000000E-16L,
  89.  7.37590000000000000000E-17L,
  90. -8.26287000000000000000E-14L,
  91.  2.50341000000000000000E-13L,
  92.  2.46508390000000000000E-11L,
  93. -5.40004410000000000000E-11L,
  94.  1.32115526000000000000E-9L,
  95. -5.99873702700000000000E-7L,
  96. -1.62427970910000000000E-5L,
  97.  2.27849553700000000000E-3L,
  98.  0.00000000000000000000E0L
  99. };
  100. #endif
  101. #if IBMPC
  102. static short pAcof[] = {
  103. 0x3761,0xf547,0x551b,0xee0b,0xbfe0, XPD
  104. 0x6a6d,0x43d6,0xc224,0xcc65,0xbfe6, XPD
  105. 0x3ffe,0x5966,0x33cb,0x8223,0x3fe9, XPD
  106. 0x9413,0x06aa,0x98c4,0xdbb2,0x3fee, XPD
  107. 0xdda4,0x9c6c,0xabe2,0xb6fc,0x3ff2, XPD
  108. 0xc6e3,0xe62d,0x86e7,0xeca6,0xbff5, XPD
  109. 0xe8c0,0xe6f2,0xad70,0xf0f6,0xbffc, XPD
  110. 0x6018,0x9d1f,0xf2ba,0x9e59,0x3ffb, XPD
  111. 0x4c98,0x8c15,0xea4a,0xde64,0x4005, XPD
  112. 0xef9e,0xc6a7,0xf74b,0xc472,0x400e, XPD
  113. };
  114. static short nodecof[] = {
  115. 0x38c8,0xf647,0x072a,0xbf64,0x3fcc, XPD
  116. 0x9b76,0xc03f,0x989c,0xc1f1,0xbfce, XPD
  117. 0x88e7,0x5905,0x4e3b,0xd9b9,0xbfd7, XPD
  118. 0xa57c,0x25f2,0xfb80,0x8450,0x3fda, XPD
  119. 0xf7b9,0x5a5a,0x1a04,0xadb2,0x3fe0, XPD
  120. 0x41a8,0xe90f,0x1783,0xef69,0xbfe2, XPD
  121. 0x8c7c,0x7312,0x2d05,0xc2ad,0xbfe8, XPD
  122. 0xe099,0x5ad6,0x1b01,0xec41,0x3fe5, XPD
  123. 0x8ff0,0x1112,0x428b,0x9c04,0x3ff1, XPD
  124. 0xbe8d,0x7207,0x9d56,0xac5a,0xbffa, XPD
  125. 0xa4e6,0x34d2,0xd051,0xc355,0x4000, XPD
  126. };
  127. static short inclcof[] = {
  128. 0x0681,0xeffb,0x9db4,0x8c0b,0x3fca, XPD
  129. 0xf1e2,0xed33,0xa0f0,0xaa13,0x3fc9, XPD
  130. 0xcd60,0x3987,0x33db,0xba10,0xbfd3, XPD
  131. 0xb113,0x604a,0xf0b7,0x8ced,0x3fd5, XPD
  132. 0x206c,0xe1ba,0xc131,0xd8d4,0x3fdb, XPD
  133. 0xb3b2,0xfa51,0x176b,0xed7f,0xbfdc, XPD
  134. 0xd03a,0x5b5d,0x04ac,0xb594,0x3fe1, XPD
  135. 0xe322,0xf2f6,0x01c7,0xa107,0xbfea, XPD
  136. 0x7f4e,0x73db,0x2422,0x8841,0xbfef, XPD
  137. 0x9c64,0xc468,0xcfd0,0x9552,0x3ff6, XPD
  138. 0x0000,0x0000,0x0000,0x0000,0x0000, XPD
  139. };
  140. #endif
  141. #if MIEEE
  142. static long pAcof[30] = {
  143. 0xbfe00000,0xee0b551b,0xf5473761,
  144. 0xbfe60000,0xcc65c224,0x43d66a6d,
  145. 0x3fe90000,0x822333cb,0x59663ffe,
  146. 0x3fee0000,0xdbb298c4,0x06aa9413,
  147. 0x3ff20000,0xb6fcabe2,0x9c6cdda4,
  148. 0xbff50000,0xeca686e7,0xe62dc6e3,
  149. 0xbffc0000,0xf0f6ad70,0xe6f2e8c0,
  150. 0x3ffb0000,0x9e59f2ba,0x9d1f6018,
  151. 0x40050000,0xde64ea4a,0x8c154c98,
  152. 0x400e0000,0xc472f74b,0xc6a7ef9e
  153. };
  154. static long nodecof[33] = {
  155. 0x3fcc0000,0xbf64072a,0xf64738c8,
  156. 0xbfce0000,0xc1f1989c,0xc03f9b76,
  157. 0xbfd70000,0xd9b94e3b,0x590588e7,
  158. 0x3fda0000,0x8450fb80,0x25f2a57c,
  159. 0x3fe00000,0xadb21a04,0x5a5af7b9,
  160. 0xbfe20000,0xef691783,0xe90f41a8,
  161. 0xbfe80000,0xc2ad2d05,0x73128c7c,
  162. 0x3fe50000,0xec411b01,0x5ad6e099,
  163. 0x3ff10000,0x9c04428b,0x11128ff0,
  164. 0xbffa0000,0xac5a9d56,0x7207be8d,
  165. 0x40000000,0xc355d051,0x34d2a4e6
  166. };
  167. static long inclcof[33] = {
  168. 0x3fca0000,0x8c0b9db4,0xeffb0681,
  169. 0x3fc90000,0xaa13a0f0,0xed33f1e2,
  170. 0xbfd30000,0xba1033db,0x3987cd60,
  171. 0x3fd50000,0x8cedf0b7,0x604ab113,
  172. 0x3fdb0000,0xd8d4c131,0xe1ba206c,
  173. 0xbfdc0000,0xed7f176b,0xfa51b3b2,
  174. 0x3fe10000,0xb59404ac,0x5b5dd03a,
  175. 0xbfea0000,0xa10701c7,0xf2f6e322,
  176. 0xbfef0000,0x88412422,0x73db7f4e,
  177. 0x3ff60000,0x9552cfd0,0xc4689c64,
  178. 0x00000000,0x00000000,0x00000000
  179. };
  180. #endif
  181.  
  182. #else /* regular double */
  183. #if UNK
  184. static DOUBLE pAcof[10] = {
  185. -8.66000000000000000000E-10,
  186. -4.75900000000000000000E-8,
  187.  2.42400000000000000000E-7,
  188.  1.30950000000000000000E-5,
  189.  1.74510000000000000000E-4,
  190. -1.80550000000000000000E-3,
  191. -2.35316000000000000000E-1,
  192.  7.73200000000000000000E-2,
  193.  1.11197100000000000000E2,
  194.  5.02909660000000000000E4
  195. };
  196. static DOUBLE nodecof[11] = {
  197.  6.64020000000000000000E-16,
  198. -2.69151000000000000000E-15,
  199. -1.54702100000000000000E-12,
  200.  7.52131300000000000000E-12,
  201.  6.31901310000000000000E-10,
  202. -3.48388152000000000000E-9,
  203. -1.81306589600000000000E-7,
  204.  2.75036225000000000000E-8,
  205.  7.43945314260000000000E-5,
  206. -4.20786043170000000000E-2,
  207.  3.05211265497500000000E0
  208. };
  209. static DOUBLE inclcof[11] = {
  210.  1.21470000000000000000E-16,
  211.  7.37590000000000000000E-17,
  212. -8.26287000000000000000E-14,
  213.  2.50341000000000000000E-13,
  214.  2.46508390000000000000E-11,
  215. -5.40004410000000000000E-11,
  216.  1.32115526000000000000E-9,
  217. -5.99873702700000000000E-7,
  218. -1.62427970910000000000E-5,
  219.  2.27849553700000000000E-3,
  220.  0.00000000000000000000E0
  221. };
  222. #endif
  223. #if DEC
  224. static short pAcof[40] = {
  225. 0130556,0005525,0015765,0043467,
  226. 0132114,0062702,0022103,0153152,
  227. 0032602,0021463,0145531,0063100,
  228. 0034133,0131230,0142006,0125224,
  229. 0035066,0176253,0161234,0066336,
  230. 0135754,0123206,0163746,0026707,
  231. 0137560,0173255,0070346,0171351,
  232. 0037236,0054762,0135235,0017540,
  233. 0041736,0062352,0045214,0012515,
  234. 0044104,0071367,0045706,0123760
  235. };
  236. static short nodecof[44] = {
  237. 0023477,0062007,0025366,0043471,
  238. 0124101,0170630,0116300,0037633,
  239. 0126331,0134516,0035531,0002611,
  240. 0027004,0050373,0100045,0171245,
  241. 0030455,0131032,0002132,0055370,
  242. 0131157,0064427,0101751,0007502,
  243. 0132502,0126455,0002563,0011214,
  244. 0031754,0040433,0000532,0153341,
  245. 0034634,0002102,0105421,0011220,
  246. 0137054,0055235,0053162,0003677,
  247. 0040503,0052720,0050464,0151245
  248. };
  249. static short inclcof[44] = {
  250. 0023014,0005635,0132357,0175406,
  251. 0022652,0011640,0170355,0031762,
  252. 0125272,0010063,0155471,0103715,
  253. 0025614,0166760,0133540,0045261,
  254. 0027330,0152301,0030741,0135040,
  255. 0127555,0077427,0065772,0050664,
  256. 0030665,0112004,0126133,0056720,
  257. 0133041,0003401,0143762,0173343,
  258. 0134210,0040444,0021163,0155577,
  259. 0036025,0051317,0150304,0064234,
  260. 0000000,0000000,0000000,0000000,
  261. };
  262. #endif
  263. #if IBMPC
  264. static short pAcof[40] = {
  265. 0xa8e7,0xa37e,0xc16a,0xbe0d,
  266. 0x7acd,0x4488,0x8cb8,0xbe69,
  267. 0x2cc8,0x796b,0x4466,0x3e90,
  268. 0xd553,0x1880,0x7653,0x3eeb,
  269. 0x8d9c,0x7c53,0xdf95,0x3f26,
  270. 0xc5b9,0xdcfc,0x94d0,0xbf5d,
  271. 0xde5d,0xae1c,0x1ed5,0xbfce,
  272. 0xa3ec,0x5753,0xcb3e,0x3fb3,
  273. 0x82aa,0x4951,0xcc9d,0x405b,
  274. 0xd4fe,0xe978,0x8e5e,0x40e8
  275. };
  276. static short nodecof[44] = {
  277. 0xc8e7,0xe55e,0xec80,0x3cc7,
  278. 0x07f3,0x1398,0x3e33,0xbce8,
  279. 0x20b1,0xc76b,0x3729,0xbd7b,
  280. 0xbe55,0x7004,0x8a1f,0x3da0,
  281. 0x4b5f,0x408b,0xb643,0x3e05,
  282. 0x21e8,0xf07d,0xed22,0xbe2d,
  283. 0x6252,0xa0ae,0x55a5,0xbe88,
  284. 0x5adc,0x602b,0x8823,0x3e5d,
  285. 0x2252,0x5162,0x8088,0x3f13,
  286. 0x40f8,0xaace,0x8b53,0xbfa5,
  287. 0x9a55,0x0a26,0x6aba,0x4008
  288. };
  289. static short inclcof[44] = {
  290. 0xff61,0xb69d,0x8173,0x3ca1,
  291. 0xa67e,0x1e1d,0x4274,0x3c95,
  292. 0x30fa,0x7b67,0x4206,0xbd37,
  293. 0x0956,0x16ec,0x9dbe,0x3d51,
  294. 0x3744,0x263c,0x1a98,0x3dbb,
  295. 0x4a36,0xed7f,0xafe2,0xbdcd,
  296. 0x6bba,0x958b,0xb280,0x3e16,
  297. 0x5edc,0x38fe,0x20e0,0xbea4,
  298. 0x7b70,0x844e,0x0824,0xbef1,
  299. 0x8d14,0xfa18,0xaa59,0x3f62,
  300. 0x0000,0x0000,0x0000,0x0000,
  301. };
  302. #endif
  303. #if MIEEE
  304. static short pAcof[40] = {
  305. 0xbe0d,0xc16a,0xa37e,0xa8e7,
  306. 0xbe69,0x8cb8,0x4488,0x7acd,
  307. 0x3e90,0x4466,0x796b,0x2cc8,
  308. 0x3eeb,0x7653,0x1880,0xd553,
  309. 0x3f26,0xdf95,0x7c53,0x8d9c,
  310. 0xbf5d,0x94d0,0xdcfc,0xc5b9,
  311. 0xbfce,0x1ed5,0xae1c,0xde5d,
  312. 0x3fb3,0xcb3e,0x5753,0xa3ec,
  313. 0x405b,0xcc9d,0x4951,0x82aa,
  314. 0x40e8,0x8e5e,0xe978,0xd4fe
  315. };
  316. static short nodecof[44] = {
  317. 0x3cc7,0xec80,0xe55e,0xc8e7,
  318. 0xbce8,0x3e33,0x1398,0x07f3,
  319. 0xbd7b,0x3729,0xc76b,0x20b1,
  320. 0x3da0,0x8a1f,0x7004,0xbe55,
  321. 0x3e05,0xb643,0x408b,0x4b5f,
  322. 0xbe2d,0xed22,0xf07d,0x21e8,
  323. 0xbe88,0x55a5,0xa0ae,0x6252,
  324. 0x3e5d,0x8823,0x602b,0x5adc,
  325. 0x3f13,0x8088,0x5162,0x2252,
  326. 0xbfa5,0x8b53,0xaace,0x40f8,
  327. 0x4008,0x6aba,0x0a26,0x9a55
  328. };
  329. static short inclcof[44] = {
  330. 0x3ca1,0x8173,0xb69d,0xff61,
  331. 0x3c95,0x4274,0x1e1d,0xa67e,
  332. 0xbd37,0x4206,0x7b67,0x30fa,
  333. 0x3d51,0x9dbe,0x16ec,0x0956,
  334. 0x3dbb,0x1a98,0x263c,0x3744,
  335. 0xbdcd,0xafe2,0xed7f,0x4a36,
  336. 0x3e16,0xb280,0x958b,0x6bba,
  337. 0xbea4,0x20e0,0x38fe,0x5edc,
  338. 0xbef1,0x0824,0x844e,0x7b70,
  339. 0x3f62,0xaa59,0xfa18,0x8d14,
  340. 0x0000,0x0000,0x0000,0x0000
  341. };
  342. #endif
  343.  
  344. #endif
  345.  
  346. extern DOUBLE J2000; /* = 2451545.0, 2000 January 1.5 */
  347. extern DOUBLE B1950;
  348. extern DOUBLE STR; /* = 4.8481368110953599359e-6 radians per arc second */
  349. extern DOUBLE Jmillenium; /* 365250.0 days */
  350. extern DOUBLE coseps, sineps; /* see epsiln.c */
  351.  
  352. /* Subroutine arguments:
  353.  *
  354.  * R = rectangular equatorial coordinate vector to be precessed.
  355.  *     The result is written back into the input vector.
  356.  * J = Julian date
  357.  * direction =
  358.  *      Precess from J to J2000: direction = 1
  359.  *      Precess from J2000 to J: direction = -1
  360.  * Note that if you want to precess from J1 to J2, you would
  361.  * first go from J1 to J2000, then call the program again
  362.  * to go from J2000 to J2.
  363.  */
  364.  
  365. precess( R, J, direction )
  366. DOUBLE R[], J;
  367. int direction;
  368. {
  369. DOUBLE A, B, T, z, pA, W;
  370. DOUBLE x[3];
  371. DOUBLE *p;
  372. int i;
  373. DOUBLE SIN(), COS(), FABS();
  374.  
  375. if( J == B1950 )
  376.     goto done;
  377.  
  378. x[0] = R[0];
  379. x[1] = R[1];
  380. x[2] = R[2];
  381.  
  382. /* if going from B1950 to date, first rotate to J2000. */
  383. if( direction == -1 )
  384.     r118_200(x);
  385.  
  386. T = (J - J2000)/Jmillenium;
  387.  
  388. /* Implementation by elementary rotations using Laskar's expansions.
  389.  * First rotate about the x axis from the initial equator
  390.  * to the ecliptic. (The input is equatorial.)
  391.  */
  392. if( direction == 1 )
  393.     epsiln( J ); /* To J2000 */
  394. else
  395.     epsiln( J2000 ); /* From J2000 */
  396. z = coseps*x[1] + sineps*x[2];
  397. x[2] = -sineps*x[1] + coseps*x[2];
  398. x[1] = z;
  399.  
  400. /* Precession in longitude
  401.  */
  402. /* T /= 10.0; */ /* thousands of years */
  403. p = (DOUBLE *)pAcof;
  404. pA = *p++;
  405. for( i=0; i<9; i++ )
  406.     pA = pA * T + *p++;
  407. pA *= STR * T;
  408.  
  409. /* Node of the moving ecliptic on the J2000 ecliptic.
  410.  */
  411. p = (DOUBLE *)nodecof;
  412. W = *p++;
  413. for( i=0; i<10; i++ )
  414.     W = W * T + *p++;
  415.  
  416. /* Rotate about z axis to the node.
  417.  */
  418. if( direction == 1 )
  419.     z = W + pA;
  420. else
  421.     z = W;
  422. B = COS(z);
  423. A = SIN(z);
  424. z = B * x[0] + A * x[1];
  425. x[1] = -A * x[0] + B * x[1];
  426. x[0] = z;
  427.  
  428. /* Rotate about new x axis by the inclination of the moving
  429.  * ecliptic on the J2000 ecliptic.
  430.  */
  431. p = (DOUBLE *)inclcof;
  432. z = *p++;
  433. for( i=0; i<10; i++ )
  434.     z = z * T + *p++;
  435. if( direction == 1 )
  436.     z = -z;
  437. B = COS(z);
  438. A = SIN(z);
  439. z = B * x[1] + A * x[2];
  440. x[2] = -A * x[1] + B * x[2];
  441. x[1] = z;
  442.  
  443. /* Rotate about new z axis back from the node.
  444.  */
  445. if( direction == 1 )
  446.     z = -W;
  447. else
  448.     z = -W - pA;
  449. B = COS(z);
  450. A = SIN(z);
  451. z = B * x[0] + A * x[1];
  452. x[1] = -A * x[0] + B * x[1];
  453. x[0] = z;
  454.  
  455. /* Rotate about x axis to final equator.
  456.  */
  457. if( direction == 1 )
  458.     epsiln( J2000 );
  459. else
  460.     epsiln( J );
  461. z = coseps * x[1] - sineps * x[2];
  462. x[2] = sineps * x[1] + coseps * x[2];
  463. x[1] = z;
  464.  
  465. /* if going from date to 1950, rotate J2000 answer to 1950. */
  466. if( direction == 1 )
  467.     r200_118(x);
  468.  
  469. for( i=0; i<3; i++ )
  470.     R[i] = x[i];
  471. done: ;
  472. }
  473.